module functions
   use iso_fortran_env, only: dp => real64
   use iso_c_binding

contains

   subroutine RHSFunction(u, rhs, t, ctx)
      real(dp), intent(in), target :: u(:)
      real(dp), intent(inout), target :: rhs(:)
      real(dp), intent(in) :: t
      type(c_ptr), intent(inout) :: ctx

      rhs(1) = -u(2)
      rhs(2) = u(1)

   end subroutine RHSFunction

   subroutine ExactSolution(t, u)
      real(dp), intent(in) :: t
      real(dp), intent(inout) :: u(:)

      u(1) = cos(t)
      u(2) = sin(t)

   end subroutine ExactSolution

   real(dp) function norm(u) result(r)
      real(dp), intent(in) :: u(:)

      r = sqrt(sum(u*u)/size(u))
   end function norm

end module functions
program scalar_eqn
   use mod_odepack
   use petsc_wrapper
   use iso_fortran_env, only: dp => real64
   use functions

   integer :: i

   type(type_ode_solver) :: ode_ts
   type(c_ptr) :: ctx
   real(dp) :: u(2), u_exact(2)
   real(dp) :: dt = 0.01_dp, t0 = 0.01_dp

   call petsc_initialize()
   call ode_ts%create("rk_tvd3")

   call ode_ts%ode_solver%set_solution(u)
   call ode_ts%ode_solver%set_time_step(dt)
   call ode_ts%ode_solver%set_rhs_function(array1d_function=RHSFunction)
   call exactSolution(t0, u)
   do i = 1, 990
      call ode_ts%ode_solver%step(u, ctx)
   end do
   call ExactSolution(t0 + (i - 1)*dt, u_exact)
   print *, u - u_exact
   print *, norm(u - u_exact)
   call petsc_finalize()
end
